# Replication file 2 

################################################
## Appendix 5. Out-of-sample test of hypotheses.
################################################

## Data preparation from CSES dataset 


###
setwd("cses4_stata") # set working directory to folder where data file is stored
library(haven)

df <- read_dta("cses4.dta") # load cses module 4 dataset 

# D1006_NAM - name of country 
df <- subset(df, D1006_NAM %in% c("France", "Germany", "Austria")) # subset of countries for comparision 


## Vote_federal_incumbent				D3006_LH_PL D3006_PR_1	

# France
# 1. D - Nathalie Arthaud = LO = L
# 2. U - Philippe Poutou = NPA = L
# 3. N - Jean-Luc Mélenchon = LF = L
# 4. H - Eva Joly = EELV = L
# 5. R - François Hollande = PS = L
# 6. S - François Bayrou = MoDem = NA
# 7. E - Nicolas Sarkozy = UPM = CR
# 8. A - Nicolas Dupont-Aignan = DLF = R
# 9. J - Marine Le Pen = FN = R
# 10. K - Jacques Cheminade = LM = NA
# 96. Cast invalid ballot (blank or null) = 
#REF
#DK 

#Germany
# 01.           Christlich Demokratische Union/Christlich Soziale Union (Union) = CR
# 02.  PARTY A  Christlich Demokratische Union Deutschlands (CDU) = CR
# 03.  PARTY E  Christlich Soziale Union in Bayern (CSU) = CR 
# 04.  PARTY B  Sozialdemokratische Partei Deutschlands (SPD) Social Democratic Party = CL                    
# 05.  PARTY C  Die Linke (Linke) = L 
# 06.  PARTY D  Buendnis90/Die Gruenen (Gruene) Alliance 90; The Greens = L
# 07.  PARTY F  Freie Demokratische Partei (FDP) = R
# 08.  PARTY G  Alternative fuer Deutschland (AfD) = R 
# 09.  PARTY H  Piratenpartei Deutschland (Piraten) Pirate Party of Germany = L 
# 10.  PARTY I  Nationaldemokratische Partei Deutschlands (NPD) = R 
# 11.           Freie Waehler (FW) = NA
# 12.           Partei Mensch Umwelt Tierschutz (Tierschutzpartei) Animal Protection Party = CL 
# 13.           Oekologisch-Demokratische Partei (ODP) Ecological Democratic German Party  = L
# 14.           Republikaner (REP) = R 
# 15.           Rentner Partei Deutschland (RENTNER) Pensioners/Greys = NA  
# 16.           Bayernpartei (BP) Bavarian Party = CR 
# 17.           Partei Bibeltreuer Christen (PBC) = NA
# 18.           Buergerrechtsbewegung Solidaritaet (BueSo) = NA
# 19.           Die Violetten (DIE VIOLETTEN) = NA
# 20.           Partei fuer Arbeit, Rechtsstaat, Tierschutz, Elitenfoerderung und basisdemokratische Initiative (Die PARTEI) = NA
# 21.           DIE GRAUEN - Graue Panther (GRAUE) = NA

#Austria
# 01.  PARTY A  Sozialdemokratische Partei Oesterreichs (SPO) Social Democratic Party of Austria = CL 
# 02.  PARTY B  Oesterreichische Volkspartei (OVP) Austrian People's Party = CR
# 03.  PARTY C  Freiheitliche Partei Oesterreichs (FPO) Freedom Party of Austria = R
# 04.  PARTY D  Gruene (GRUENE) Green Party = CL
# 05.  PARTY G  Team Stronach fuer Oesterreich (TEAM STRONACH) = R 
# 06.  PARTY F  NEOS - Das Neue Oesterreich und Liberales Forum (NEOS) - The New Austria = CR
# 07.  PARTY E  Buendnis Zukunft Oesterreich (BZO) Alliance for the Future of Austria = R
# 08.  PARTY H  Kommunistische Partei Oesterreichs (KPO) - Communist Party of Austria = L
# 09.  PARTY I  Piratenpartei Osterreichs (PIRATEN) - SEE ELECTION STUDY NOTE Pirate Party = ?
# 10.           Buergerforum Oesterreich (FRITZ) Citizens' Forum Austria = L 
# 11.           Liberales Forum (LiF) Liberal Forum
# 12.           Die Freiheitlichen in Kaernten (FPK) Freedom Party in Carinthia


## Evolution_economy					D3003_1
## Patrimony_high						D3027_2		D3027_3
## Patrimony_low						D3027_1		D3027_4
## Income_government_action				D3004	
## Trust_fedparl	
## trust_index	
## Europeanisation	
## immigration_index	
## Party_close							D3018_1
## Leftright_self						D3014
## Interest_politics_general	
## Satisfaction_fedgov	
## Age									D2001_Y year of birth 
## Sex									D2002
## Marital								D2004
## Education							D2003
## Employment							D2010
## Religion_frequency					D2024


## Subset of variables used in the analysis
df1 <- df[, c("D1004", "D3006_LH_PL", "D3006_PR_1", "D3003_1", "D3027_2", "D3027_3", "D3027_1", "D3027_4", 
			  "D3004", "D3018_1", "D3014", "D2001_Y", "D2002", "D2004", "D2003", "D2010", "D2024", "D5007", "D5008")]

## Recodeing variables 

#    D5007        >>> M02a    PARTY OF THE PRESIDENT BEFORE
#    D5008        >>> M02b    PARTY OF THE PRIME MINISTER BEFORE 
#	 D1004 		  >>> ID VARIABLE - ELECTION STUDY (ALPHABETIC POLITY)
df1$country <- as.numeric(as.factor(df1$D1004))

df1$Vote_federal_incumbent_Germany <- ifelse(df1$country == 2, df1$D3006_LH_PL, NA)
df1$Vote_federal_incumbent_Austria <- ifelse(df1$country == 1, df1$D3006_LH_PL, NA)
df1$Vote_federal_incumbent_France <- df1$D3006_PR_1

df1$Evolution_economy <- df1$D3003_1
	df1$Evolution_economy <- ifelse(df1$Evolution_economy == 5, 1,
                                       ifelse(df1$Evolution_economy == 3, 2,
                                              ifelse(df1$Evolution_economy == 1, 3, NA)))

df1$Patrimony_business <- df1$D3027_2
	df1$Patrimony_business <- ifelse(df1$Patrimony_business == 5, 0, 
									ifelse(df1$Patrimony_business > 5, NA, df1$Patrimony_business))

df1$Patrimony_stocks <- df1$D3027_3		
	df1$Patrimony_stocks <- ifelse(df1$Patrimony_stocks == 5, 0, 
									ifelse(df1$Patrimony_stocks > 5, NA, df1$Patrimony_stocks))

df1$Patrimony_house <- df1$D3027_1	
	df1$Patrimony_house <- ifelse(df1$Patrimony_house == 5, 0, 
									ifelse(df1$Patrimony_house > 5, NA, df1$Patrimony_house))	

df1$Patrimony_savings <- df1$D3027_4	
	df1$Patrimony_savings <- ifelse(df1$Patrimony_savings == 5, 0, 
									ifelse(df1$Patrimony_savings > 5, NA, df1$Patrimony_savings))	

df1$Income_government_action <- df1$D3004
	df1$Income_government_action <- ifelse(df1$Income_government_action > 5, NA, df1$Income_government_action)

df1$Party_close <- df1$D3018_1
	df1$Party_close <- ifelse(df1$Party_close > 5, NA,
							ifelse(df1$Party_close == 5, 0, df1$Party_close))

df1$Leftright_self <- df1$D3014
	df1$Leftright_self <- ifelse(df1$Leftright_self > 10, NA, df1$Leftright_self)

df1$Age <- df1$D2001_Y
	df1$Age <- ifelse(df1$Age > 5000, NA, df1$Age)

df1$Sex <- df1$D2002
	df1$Sex <- ifelse(df1$Sex == 2, 0, df1$Sex)

df1$Marital <- df1$D2004
	df1$Marital <- ifelse(df1$Marital > 5, NA, 
						ifelse(df1$Marital > 1, 0, df1$Marital))

df1$Education <- df1$D2003
	df1$Education <- ifelse(df1$Education > 10, NA, df1$Education)

df1$Employment <- df1$D2010
	df1$Employment <- ifelse(df1$Employment > 15, NA, df1$Employment)

df1$Religion_frequency <- df1$D2024
	df1$Religion_frequency <- ifelse(df1$Religion_frequency > 6, NA, df1$Religion_frequency)

df1$Austria_LR <- ifelse(df1$Vote_federal_incumbent_Austria > 10, NA, 
					(ifelse(df1$Vote_federal_incumbent_Austria == 1, 0, 
						(ifelse(df1$Vote_federal_incumbent_Austria == 2, 1,
							(ifelse(df1$Vote_federal_incumbent_Austria == 3, 1,
								(ifelse(df1$Vote_federal_incumbent_Austria == 4, 0,
									(ifelse(df1$Vote_federal_incumbent_Austria == 5, 1,
										(ifelse(df1$Vote_federal_incumbent_Austria == 6, 1,
											(ifelse(df1$Vote_federal_incumbent_Austria == 7, 1,
												(ifelse(df1$Vote_federal_incumbent_Austria == 8, 0,
													(ifelse(df1$Vote_federal_incumbent_Austria == 9, 0,
														(ifelse(df1$Vote_federal_incumbent_Austria == 10, NA, df1$Vote_federal_incumbent_Austria)))))))))))))))))))))


df1$France_LR <- ifelse(df1$Vote_federal_incumbent_France > 9, NA,
					(ifelse(df1$Vote_federal_incumbent_France == 6, NA,
					 	(ifelse(df1$Vote_federal_incumbent_France < 6, 0,
					 		(ifelse(df1$Vote_federal_incumbent_France > 6, 1, df1$Vote_federal_incumbent_France)))))))

df1$Germany_LR <- ifelse(df1$Vote_federal_incumbent_Germany > 21, NA,
				 	(ifelse(df1$Vote_federal_incumbent_Germany == 1, 1,
				 		(ifelse(df1$Vote_federal_incumbent_Germany == 4, 0,
				 			(ifelse(df1$Vote_federal_incumbent_Germany == 5, 0,
				 				(ifelse(df1$Vote_federal_incumbent_Germany == 6, 0,
				 					(ifelse(df1$Vote_federal_incumbent_Germany == 7, 1,
				 						(ifelse(df1$Vote_federal_incumbent_Germany == 8, 1,
				 							(ifelse(df1$Vote_federal_incumbent_Germany == 9, 0,
				 								(ifelse(df1$Vote_federal_incumbent_Germany == 10, 1,
				 									(ifelse(df1$Vote_federal_incumbent_Germany == 11, 1,
				 										(ifelse(df1$Vote_federal_incumbent_Germany == 12, 0,
				 											(ifelse(df1$Vote_federal_incumbent_Germany == 13, 0,
				 												(ifelse(df1$Vote_federal_incumbent_Germany == 14, 1,
				 													(ifelse(df1$Vote_federal_incumbent_Germany == 15, NA,
				 														(ifelse(df1$Vote_federal_incumbent_Germany == 17, NA,
				 															(ifelse(df1$Vote_federal_incumbent_Germany == 20, NA, df1$Vote_federal_incumbent_Germany)))))))))))))))))))))))))))))))


df1$LR_parl <- ifelse(is.na(df1$Germany_LR), df1$Austria_LR, df1$Germany_LR) 

df1$Patrimony_high <- df1$Patrimony_business + df1$Patrimony_stocks 
df1$Patrimony_low <- df1$Patrimony_house + df1$Patrimony_savings

## End of recoding variables. 

####################################################################################################
## eTable 5a. Logistic regression results for out-of-sample test of Multidimensional economic voting
####################################################################################################

## Germany
modGerLR1 <- glm(Germany_LR ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + Party_close + Leftright_self + Age + Sex + Marital + Education + Employment + Religion_frequency, data = df1, family = binomial()) 

## France
modFraLR1 <- glm(France_LR ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + Party_close + Leftright_self + Age + Sex + Education + Employment + Religion_frequency, data = df1, family = binomial()) 

# Note Marital status in France was missing and that is why it is excluded. No observations. 

## Austria
modAusLR1 <- glm(Austria_LR ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + Party_close + Leftright_self + Age + Sex + Marital + Education + Religion_frequency, data = df1, family = binomial())

# Note Employment in Austria was differently asked and misses half of the observations, that is why it is excluded. 

# Pooled sample 
modParlLR1 <- glm(LR_parl ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + Party_close + Leftright_self + Age + Sex + Marital + Education + Religion_frequency + as.factor(country), data = df1, family = binomial())

###################################################################################################
## eTable 5b. Linear Probability Models with White Corrected Standard Errors for Heteroscedasticity
###################################################################################################

## Germany LPM
modGerLR1_lpm <- lm(Germany_LR ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + Party_close + Leftright_self + Age + Sex + Marital + Education + Employment + Religion_frequency, data = df1) 

## France LPM
modFraLR1_lpm <- lm(France_LR ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + Party_close + Leftright_self + Age + Sex + Education + Employment + Religion_frequency, data = df1) 

# Note Marital status in France was missing and that is why it is excluded. No observations. 

## Austria LPM
modAusLR1_lpm <- lm(Austria_LR ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + Party_close + Leftright_self + Age + Sex + Marital + Education + Religion_frequency, data = df1)

# Note Employment in Austria was differently asked and misses half of the observations, that is why it is excluded. 

# Pooled sample 
modParlLR1_lpm <- lm(LR_parl ~ Evolution_economy + Patrimony_high + Patrimony_low + Income_government_action + Party_close + Leftright_self + Age + Sex + Marital + Education + Religion_frequency + as.factor(country), data = df1)

library(sandwich)
library(lmtest)

## Corrected standard errors
vv1 <- vcovHC(modParlLR1_lpm, type="HC1")
coeftest(modParlLR1_lpm, vcov = vv1)

vv1a <- vcovHC(modGerLR1_lpm, type="HC1")
coeftest(modGerLR1_lpm, vcov = vv1a)

vv1b <- vcovHC(modAusLR1_lpm, type="HC1")
coeftest(modAusLR1_lpm, vcov = vv1b)

vv2 <- vcovHC(modFraLR1_lpm, type="HC1")
coeftest(modFraLR1_lpm, vcov = vv2)


#################################################################################################################
## eFigure 5a. Multidimensional Voting in Germany and Auestia compared to Flanders regional and federal elecitons
#################################################################################################################
library(ggplot2)

##############################################################
## Pooled model (parliamentary elections) top plot - left side
##############################################################

library(margins)
marginal_effects <- margins(modParlLR1)
marginal_effects <- margins(modParlLR1, variables = c("Evolution_economy", "Patrimony_high", "Patrimony_low", "Income_government_action"))
summary(marginal_effects)

#                   factor    AME     SE      z      p   lower  upper
#        Evolution_economy 0.0298 0.0168 1.7805 0.0750 -0.0030 0.0627
# Income_government_action 0.0419 0.0094 4.4287 0.0000  0.0233 0.0604
#           Patrimony_high 0.0506 0.0192 2.6410 0.0083  0.0130 0.0881
#            Patrimony_low 0.0170 0.0192 0.8867 0.3752 -0.0206 0.0547

# Create a data frame with your summary data
summary_data <- data.frame(
  Variable = factor(c("Economic evaluation", "Patrimony (high)", "Patrimony (low)", "Positional Economics"),
                    levels = c("Economic evaluation", "Patrimony (high)", "Patrimony (low)", "Positional Economics")),
  AME = c(0.0298, 0.0506, 0.0170, 0.0419),
  SE = c(0.0168, 0.0192, 0.0192, 0.0094),
  p_value = c(0.0750, 0.0083, 0.3752, 0.0000)  
)

# Function to add stars for significance levels next to AME
add_stars <- function(AME, p_value) {
  stars <- ifelse(p_value < 0.001, "***",
           ifelse(p_value < 0.01, "**",
           ifelse(p_value < 0.05, "*", 
           ifelse(p_value < 0.1, "+", ""))))
  return(paste0(sprintf("%.3f", AME), stars))
}

# Add AME with significance stars to the data frame
summary_data$AME_with_Stars <- mapply(add_stars, summary_data$AME, summary_data$p_value)

# Plot the AME with error bars for SE and display rounded AME values with significance stars
modParlLR1_ame <- ggplot(data = summary_data, aes(x = Variable, y = AME, ymin = AME - SE, ymax = AME + SE)) +
  geom_bar(stat = "identity", fill = "gray", alpha = 0.5) +
  geom_errorbar(width = 0.3) + 
  geom_text(aes(label = AME_with_Stars, hjust = 1.05, vjust = -0.5), size = 8, family = "Times New Roman", fontface = "bold") +
  labs(title = "",
       x = "",
       y = "") +
  theme_minimal() +
  theme(  axis.title.y = element_text(size = 22, family = "Times New Roman", face = "bold"),
          axis.text.y = element_text(size = 22, family = "Times New Roman", face = "bold"),
          axis.title.x = element_text(size = 22, family = "Times New Roman", face = "bold"),
          axis.text.x = element_text(size = 22, family = "Times New Roman", face = "bold")) +
  scale_y_continuous(limits = c(-0.03, 0.10))

ggsave("modParlLR1_ame.png", modParlLR1_ame, width = 13, height = 8, units = "in", dpi = 300)

############################################################
## Germany (parliamentary elections) middle plot - left side
############################################################

marginal_effects <- margins(modGerLR1)
marginal_effects <- margins(modGerLR1, variables = c("Evolution_economy", "Patrimony_high", "Patrimony_low", "Income_government_action"))
summary(marginal_effects)

#                   factor    AME     SE      z      p   lower  upper
#        Evolution_economy 0.0606 0.0205 2.9589 0.0031  0.0204 0.1007
# Income_government_action 0.0312 0.0121 2.5719 0.0101  0.0074 0.0550
#           Patrimony_high 0.0524 0.0281 1.8677 0.0618 -0.0026 0.1075
#            Patrimony_low 0.0244 0.0234 1.0415 0.2976 -0.0215 0.0703

# Create a data frame with your summary data
summary_data <- data.frame(
  Variable = factor(c("Economic evaluation", "Patrimony (high)", "Patrimony (low)", "Positional Economics"),
                    levels = c("Economic evaluation", "Patrimony (high)", "Patrimony (low)", "Positional Economics")),
  AME = c(0.0606, 0.0524, 0.0244, 0.0312),
  SE = c(0.0205, 0.0281, 0.0234, 0.0121),
  p_value = c(0.0031, 0.0618, 0.2976, 0.0101)  
)

# Function to add stars for significance levels next to AME
add_stars <- function(AME, p_value) {
  stars <- ifelse(p_value < 0.001, "***",
           ifelse(p_value < 0.01, "**",
           ifelse(p_value < 0.05, "*", 
           ifelse(p_value < 0.1, "+", ""))))
  return(paste0(sprintf("%.3f", AME), stars))
}

# Add AME with significance stars to the data frame
summary_data$AME_with_Stars <- mapply(add_stars, summary_data$AME, summary_data$p_value)

# Plot the AME with error bars for SE and display rounded AME values with significance stars
modGerLR1_ame <- ggplot(data = summary_data, aes(x = Variable, y = AME, ymin = AME - SE, ymax = AME + SE)) +
  geom_bar(stat = "identity", fill = "gray", alpha = 0.5) +
  geom_errorbar(width = 0.3) + 
  geom_text(aes(label = AME_with_Stars, hjust = 1.05, vjust = -0.5), size = 8, family = "Times New Roman", fontface = "bold") +
  labs(title = "",
       x = "",
       y = "") +
  theme_minimal() +
  theme(  axis.title.y = element_text(size = 22, family = "Times New Roman", face = "bold"),
          axis.text.y = element_text(size = 22, family = "Times New Roman", face = "bold"),
          axis.title.x = element_text(size = 22, family = "Times New Roman", face = "bold"),
          axis.text.x = element_text(size = 22, family = "Times New Roman", face = "bold")) +
  scale_y_continuous(limits = c(-0.03, 0.10))

ggsave("modGerLR1_ame.png", modGerLR1_ame, width = 13, height = 8, units = "in", dpi = 300) 


############################################################
## Austria (parliamentary elections) bottom plot - left side
############################################################

marginal_effects <- margins(modAusLR1)
marginal_effects <- margins(modAusLR1, variables = c("Evolution_economy", "Patrimony_high", "Patrimony_low", "Income_government_action"))
summary(marginal_effects)

#                   factor     AME     SE       z      p   lower  upper
#        Evolution_economy -0.0250 0.0286 -0.8741 0.3821 -0.0811 0.0311
# Income_government_action  0.0496 0.0150  3.3186 0.0009  0.0203 0.0789
#           Patrimony_high  0.0501 0.0269  1.8668 0.0619 -0.0025 0.1028
#            Patrimony_low -0.0056 0.0331 -0.1690 0.8658 -0.0704 0.0592
#


# Create a data frame with your summary data
summary_data <- data.frame(
  Variable = factor(c("Economic evaluation", "Patrimony (high)", "Patrimony (low)", "Positional Economics"),
                    levels = c("Economic evaluation", "Patrimony (high)", "Patrimony (low)", "Positional Economics")),
  AME = c(-0.0250, 0.0501, -0.0056, 0.0496),
  SE = c(0.0286, 0.0269, 0.0331, 0.0150),
  p_value = c(0.3821, 0.0619, 0.8658, 0.0009)  
)

# Function to add stars for significance levels next to AME
add_stars <- function(AME, p_value) {
  stars <- ifelse(p_value < 0.001, "***",
           ifelse(p_value < 0.01, "**",
           ifelse(p_value < 0.05, "*", 
           ifelse(p_value < 0.1, "+", ""))))
  return(paste0(sprintf("%.3f", AME), stars))
}

# Add AME with significance stars to the data frame
summary_data$AME_with_Stars <- mapply(add_stars, summary_data$AME, summary_data$p_value)

# Plot the AME with error bars for SE and display rounded AME values with significance stars
modAusLR1_ame <- ggplot(data = summary_data, aes(x = Variable, y = AME, ymin = AME - SE, ymax = AME + SE)) +
  geom_bar(stat = "identity", fill = "gray", alpha = 0.5) +
  geom_errorbar(width = 0.3) + 
  geom_text(aes(label = AME_with_Stars, hjust = 1.05, vjust = -0.5), size = 8, family = "Times New Roman", fontface = "bold") +
  labs(title = "",
       x = "",
       y = "") +
  theme_minimal() +
  theme(  axis.title.y = element_text(size = 22, family = "Times New Roman", face = "bold"),
          axis.text.y = element_text(size = 22, family = "Times New Roman", face = "bold"),
          axis.title.x = element_text(size = 22, family = "Times New Roman", face = "bold"),
          axis.text.x = element_text(size = 22, family = "Times New Roman", face = "bold")) +
  scale_y_continuous(limits = c(-0.06, 0.10))

ggsave("modAusLR1_ame.png", modAusLR1_ame, width = 13, height = 8, units = "in", dpi = 300)

## Please note: the plots on the right side are available in "replication file" uder Main Analysis, Figure 2

####################################################################################################
## eFigure 5b. Multidimensional Voting in France compared to Wallonia regional and federal elecitons
####################################################################################################

##############################################
## France (Presidential elections) - left side
##############################################

marginal_effects <- margins(modFraLR1)
marginal_effects <- margins(modFraLR1, variables = c("Evolution_economy", "Patrimony_high", "Patrimony_low", "Income_government_action"))
summary(marginal_effects)

#                   factor     AME     SE       z      p   lower  upper
#        Evolution_economy -0.0094 0.0182 -0.5164 0.6056 -0.0451 0.0263
# Income_government_action  0.0444 0.0101  4.4069 0.0000  0.0246 0.0641
#           Patrimony_high  0.0148 0.0177  0.8328 0.4050 -0.0200 0.0496
#            Patrimony_low  0.0013 0.0126  0.1003 0.9201 -0.0234 0.0259

# Create a data frame with your summary data
summary_data <- data.frame(
  Variable = factor(c("Economic evaluation", "Patrimony (high)", "Patrimony (low)", "Positional Economics"),
                    levels = c("Economic evaluation", "Patrimony (high)", "Patrimony (low)", "Positional Economics")),
  AME = c(-0.0094, 0.0148, 0.0013, 0.0444),
  SE = c(0.0182, 0.0177, 0.0126, 0.0101),
  p_value = c(0.6056, 0.4050, 0.9201, 0.0000)  
)

# Function to add stars for significance levels next to AME
add_stars <- function(AME, p_value) {
  stars <- ifelse(p_value < 0.001, "***",
           ifelse(p_value < 0.01, "**",
           ifelse(p_value < 0.05, "*", 
           ifelse(p_value < 0.1, "+", ""))))
  return(paste0(sprintf("%.3f", AME), stars))
}

# Add AME with significance stars to the data frame
summary_data$AME_with_Stars <- mapply(add_stars, summary_data$AME, summary_data$p_value)

# Plot the AME with error bars for SE and display rounded AME values with significance stars
modFraLR1_ame <- ggplot(data = summary_data, aes(x = Variable, y = AME, ymin = AME - SE, ymax = AME + SE)) +
  geom_bar(stat = "identity", fill = "gray", alpha = 0.5) +
  geom_errorbar(width = 0.3) + 
  geom_text(aes(label = AME_with_Stars, hjust = 1.05, vjust = -0.5), size = 8, family = "Times New Roman", fontface = "bold") +
  labs(title = "",
       x = "",
       y = "") +
  theme_minimal() +
  theme(  axis.title.y = element_text(size = 22, family = "Times New Roman", face = "bold"),
          axis.text.y = element_text(size = 22, family = "Times New Roman", face = "bold"),
          axis.title.x = element_text(size = 22, family = "Times New Roman", face = "bold"),
          axis.text.x = element_text(size = 22, family = "Times New Roman", face = "bold")) +
  scale_y_continuous(limits = c(-0.03, 0.10))

ggsave("modFraLR1_ame.png", modFraLR1_ame, width = 13, height = 8, units = "in", dpi = 300)

## Please note: the plots on the right side are available in "replication file" uder Main Analysis, Figure 3

## End of Appendix 5. 

### End of Replication file 2.